import holoviews as hv
hv.extension('bokeh')
hv.opts.defaults(hv.opts.Curve(width=500),
hv.opts.Points(width=500),
hv.opts.Image(width=500, colorbar=True, cmap='Viridis'))
import numpy as np
import scipy.signal
import scipy.linalg
11. Estimadores adaptivos parte I¶
Hasta ahora hemos estudiando sistemas lineales donde:
sus coeficientes quedan fijos luego del diseño y son constantes en el tiempo
hacen supuestos sobre los estadísticos de la señal/ruido
¿Qué podemos hacer si
no podemos hacer supuestos sobre los estadísticos?
los estadísticos de la señal/ruido cambian en el tiempo?
estamos en un escenario donde los datos llegan continuamente (data streaming)?
Cuando los estadísticos cambian en el tiempo decimos que la señal es no estacionaria.
En estos casos necesitamos un estimador de tipo adaptivo, es decir sistemas y filtros cuyos coeficientes se pueden adaptar a medida que llegan nuevos datos. Estos estimadores se diseñan de acuerdo a un método de optimización que es online
La siguiente figura muestra algunos ejemplos de aplicaciones de sistemas adaptivos
El método de optimización online que utilizaremos principalmente en este curso es el gradiente descendente estocástico. Revisemos a continuación los fundamentos.
11.1. Gradiente descendente¶
Sea un vector de pesos \(w\) de largo \(L+1\) que guarda los coeficientes de un estimador
Sea ahora una función de costo que mapea el vector de pesos a un número real
La función de costo debe ser tal que a menor \(J\) menor sea el error del estimador
Para entrenar un estimador o filtro adaptivo se tienen los siguientes pasos conceptuales
Partimos de una solución inicial \(w_0\)
Modificamos iterativamente \(w\) tal que \(J(w_{t+1}) < J(w_t)\)
Nos detenemos al cumplir un cierto criterio
Para modificar iterativamete y eficientemente los pesos utilizaremos la regla del gradiente descendente (GD)
donde \(\mu\) se conoce como tasa de aprendizaje o “paso”
Imaginemos \(J\) como una superficie de \(L+1\) dimensiones
En cada punto el gradiente negativo de \(J\) nos indica hacia donde está el descenso más abrupto
La tasa \(\mu\) nos da el largo del salto entre \(w_t\) y \(w_{t+1}\)
Observando la expansión de Taylor de primer orden de \(J\) en \(w_{t}\)
es decir que usando la regla GD con \(\mu>0\) y asumiendo que \(J\) es convexo entonces se cumple que \(J\) siempre decrece monotónicamente.
La siguiente gráficas interactivas muestran una superficie de costo no convexa para un parámetro unidimensional. Cada punto representa una solución que parte desde una posición inicial distinta. Las flechas corresponden a la derivada multiplicada por la tasa de aprendizaje.
Estudie la evolución de las tres soluciones en cada caso. En primer lugar se utiliza \(\mu=0.05\)
J = lambda w : (w-1)**2 + 0.2*np.sin(2*np.pi*w) # Función de costo
gradJ = lambda w : 2*(w-1) + 0.2*2*np.pi*np.cos(2*np.pi*w) # Gradiente
mu = 0.05 # Tasa de aprendizaje
iteraciones = 15
wt = np.zeros(shape=(iteraciones, 3))
wt[0, :] = np.array([0.05, 0.4, 1.9]) # Solución inicial
w_plot = np.linspace(0, 2, num=100)
for k in range(1, iteraciones):
wt[k, :] = wt[k-1, :] - mu*gradJ(wt[k-1, :])
loss_surface = hv.Curve((w_plot, J(w_plot)), 'w', 'J')
hMap = hv.HoloMap(kdims='Iteración')
for k in range(iteraciones):
dots = hv.Points((wt[k, :], J(wt[k, :]))).opts(size=10, color='k')
mag = mu*gradJ(wt[k, :])
angle = np.pi/2 - np.sign(-mag)*np.pi/2
mag = np.abs(mag)
arrows = hv.VectorField((wt[k, :], J(wt[k, :]), angle, mag)).opts(pivot='tail',
magnitude=hv.dim('Magnitude'),
rescale_lengths=False)
hMap[k] = dots * arrows
loss_surface * hMap
Advertencia
Dependiendo de donde partimos la solución final es distinta. El gradiente descedente puede quedarse “atorado” en un mínimo local o en un punto silla
Ahora observe como evolucionan las tres soluciones con \(\mu=0.5\), es decir 10 veces más grande que el caso anterior
J = lambda w : (w-1)**2 + 0.2*np.sin(2*np.pi*w) # Función de costo
gradJ = lambda w : 2*(w-1) + 0.2*2*np.pi*np.cos(2*np.pi*w) # Gradiente
mu = 0.5 # Tasa de aprendizaje
iteraciones = 15
wt = np.zeros(shape=(iteraciones, 3))
wt[0, :] = np.array([0.05, 0.4, 1.9]) # Solución inicial
w_plot = np.linspace(0, 2, num=100)
for k in range(1, iteraciones):
wt[k, :] = wt[k-1, :] - mu*gradJ(wt[k-1, :])
loss_surface = hv.Curve((w_plot, J(w_plot)), 'w', 'J')
hMap = hv.HoloMap(kdims='Iteración')
for k in range(iteraciones):
dots = hv.Points((wt[k, :], J(wt[k, :]))).opts(size=10, color='k')
mag = mu*gradJ(wt[k, :])
angle = np.pi/2 - np.sign(-mag)*np.pi/2
mag = np.abs(mag)
arrows = hv.VectorField((wt[k, :], J(wt[k, :]), angle, mag)).opts(pivot='tail',
magnitude=hv.dim('Magnitude'),
rescale_lengths=False)
hMap[k] = dots * arrows
loss_surface * hMap
Advertencia
Si la tasa de aprendizaje es muy alta, los pasos son muy largos y podríamos no converger a un punto estacionario
Los ejemplos anteriores nos han mostrado algunas de las limitaciones del algoritmo de gradiente descendente. Es importante tenerlas en cuenta cuando lo utilicemos en nuestras aplicaciones
11.2. Gradiente descendente en el filtro de Wiener¶
Para el filtro de Wiener teníamos que
por ende
y finalmente
En este caso la condición para una convergencia estable es
donde \(\lambda_{\text{max}}\) es valor propio más grande de \(R_{uu}\)
(La prueba de esto puede encontrarse en Haykin, “Adaptive filter theory”, Sección 4.3)
11.3. Gradiente descendente estocástico (SGD)¶
El filtro de Wiener es óptimo pero no adaptivo:
Requiere de \(N\) muestras de \(u\) y \(d\) para estimar \(R_{ud}\) y \(R_{uu}\)
Los pesos se adaptan luego de haber presentado las \(N\) muestras: Es una estrategia de tipo batch
Asume que la señal es estacionaria
Si nuestros son no estacionarios significa que debemos adaptar el filtro a medida que nuevas muestras son observadas . Para lograr esto podemos usar la versión estocástica del GD: SGD
En SGD:
los pesos se adaptan luego de haber presentado una sola muestra o un conjunto pequeño de muestras (mini-batch)
no hay garantía de llegar al óptimo en un problema convexo, pero es más eficiente computacionalmente que GD
El siguiente esquema muestra una comparación entre la trayectoria de \(w\) cuando se usa GD (negro) y SGD (rojo). En general la trayectoria de SGD será más ruidosa y también podría requerir más pasos, pero cada paso es mucho más económico
11.4. Algoritmo Least Mean Square (LMS)¶
Podemos extender el filtro de Wiener al caso no-estacionario usando SGD, el resultado es un algoritmo simple que además es robusto: El algoritmo LMS
Fue fue inventado en 1960 por Bernard Widrow y Ted Hoff
A diferencia del filtro de Wiener no se requiere conocimiento estadístico del proceso. Tampoco se requiere calcular e invertir la matriz de correlación
El algoritmo LMS se ajusta o entrena de manera recursiva y online
Consideremos la función de costo estocástica para la arquitectura FIR que utilizamos para el filtro de Wiener
donde definimos \(\textbf{u}_n = [u_n, u_{n-1}, \ldots, u_{n-L}]\).
Nota
A diferencia del filtro de Wiener no aplicamos el valor esperado al error cuadrático. Se usa el error cuadrático instantaneo
Para continuar calculamos el gradiente en función del peso \(w_{n, k}\)
Luego, usando la regla SGD llegamos a
y que en forma matricial es
que se conoce como la regla de Widrow-Hoff
Importante
El algoritmo LMS estima el error instantaneo y actualiza los pesos recursivamente
La complejidad de este algoritmo es \(L+1\).
11.4.1. Convergencia del algoritmo LMS (Haykin 6.5)¶
El algoritmo LMS tiende en la media al valor óptimo
para \(n\to \infty\)
Además convergence en la media cuadrada: La varianza de \(\textbf{w}_n - \textbf{w}^*\) tiene al valor mínimo de \(J\) para \(n\to \infty\)
Esto se cumple si
donde \(R_{uu} = \mathbb{E}[\textbf{u}_n \textbf{u}_n^T ]\) es la matriz de autocorrelación y \(\text{Tr}[]\) el operador que calcula la traza de una matriz
11.4.2. Algoritmo Normalized LMS (NLMS)¶
Tenemos la siguiente regla iterativa
que se puede interpretar graficamente como
(donde \(\textbf{x}(k)\) y \(\textbf{w}(k)\) corresponden a \(\textbf{u}_n\) y \(\textbf{w}_n\) en nuestra notación, respectivamente)
Nota
Los cambios en el vector de peso \(\Delta \textbf{w}_n\) son paralelos a \(\textbf{u}_{n}\). Además estos cambios podrían estar dominados por
El algoritmo Normalized LMS (NLMS) corrige este problema ponderando por la varianza de \(\textbf{u}_{n}\)
donde la constante \(\delta\) es un valor pequeño que se usa para evitar divisiones por cero. En lo que sigue usaremos NLMS para revisar algunas aplicaciones
11.5. Implementación del filtro NLMS en Python¶
Podemos implementar las ecuaciones del filtro NLMS como se muestra a continuación
class Filtro_NLMS:
def __init__(self, L, mu, delta=1e-6, winit=None):
self.L = L
self.w = np.zeros(shape=(L+1, ))
self.mu = mu
self.delta = delta
def update(self, un, dn):
# Asumiendo que un = [u[n], u[n-1], ..., u[n-L]]
unorm = np.dot(un, un) + self.delta
yn = np.dot(self.w, un)
self.w += 2*self.mu*(dn - yn)*(un/unorm)
return yn
El filtro recibe como entrada el orden \(L\) y la tasa de aprendizaje \(\mu\)
Se asume un vector cero para los pesos iniciales, pero también en la práctica podemos partir de una solución anterior si esta existiera
Para actualizar el vector de pesos es necesario entregar el vector \(\textbf{u}_n \in \mathbb{R}^{L+1}\) y la salida deseada \(d_n \in \mathbb{R}\). La función
updateretorna la salida predicha por el filtro \(y_n = w_n^T \textbf{u}_n \)
A continuación probaremos este filtro con una aplicación conocida como Adaptive line enhancement (ALE). ALE se refiere a un sistema adaptivo para eliminar ruido blanco aditivo de una señal. El sistema aprende un filtro pasabanda en torno a la frecuencia de interés
En ALE usamos como señal deseada
El valor predicho por el filtro será la señal \(u\) pero libre de ruido blanco. Esto se debe a que el ruido blanco no tiene correlación y por ende el filtro adaptivo no lo puede predecir
# Digamos que u = s + n
# El objetivo es limpiar u para obtener s
# s es una señal determínista y n es ruido blanco
Fs, f0 = 100, 5
t = np.arange(0, 4, 1/Fs)
s = np.sin(2.0*np.pi*t*f0)
n = 0.5*np.random.randn(len(t))
s[t>2.0] += 5 # Simulemos un cambio abrupto en la media de la señal
#s += s*(0.5 + 0.5*np.cos(2.0*np.pi*t/2)) # Tremolo (AM)
u = s + n
A diferencia de un filtro estático (como el filtro de Wiener) es posible filtrar incluso ante cambios bruscos en la señal.
Estudie como cambia el resultado del filtro con distintos valores de \(\mu\)
L = 20
u_preds = {}
for mu in np.logspace(-2, 0, num=10):
myfilter = Filtro_NLMS(L=L, mu=mu)
u_preds[mu] = np.zeros(shape=(len(u),))
for k in range(L+1, len(u)):
u_preds[mu][k] = myfilter.update(u[k-L-1:k][::-1], u[k])
hMap = hv.HoloMap(kdims='mu')
for mu, u_pred in u_preds.items():
s1 = hv.Curve((t, s), 'Tiempo', 'Señal', label='Limpia')
s2 = hv.Scatter((t, u), 'Tiempo', 'Señal', label='Contaminada')
s3 = hv.Curve((t, u_pred), 'Tiempo', 'Señal', label='Filtrada')
hMap[mu] = hv.Overlay([s1, s2, s3]).opts(hv.opts.Overlay(legend_position='top'),
hv.opts.Curve(ylim=(-5, 10), height=350))
hMap
Importante
La tasa de aprendizaje \(\mu\) controla la velocidad de adaptación. Pero una tasa demasiado grande provoca que el filtro sea inestable. En general el valor óptimo de \(\mu\) depende del problema y del valor de \(L\)
La siguiente figura muestra la respuesta en frecuencia del filtro en función del tiempo para \(\mu=0.02\)
Observe como a medida que se adapta el filtro se concentra en la frecuencia fundamental de la señal, que en este caso es 5 Hz
L = 20
u_preds = {}
myfilter = Filtro_NLMS(L=L, mu=0.02)
H_history = np.zeros(shape=(512, len(u)))
for k in range(L+1, len(u)):
myfilter.update(u[k-L-1:k][::-1], u[k])
fk, Hk = scipy.signal.freqz(b=myfilter.w, a=1, fs=Fs)
H_history[:, k] = np.abs(Hk)
hv.Image((t, fk, H_history), kdims=['Tiempo [s]', 'Frecuencia [Hz]']).opts(cmap='Blues')
11.6. Comparación entre Filtro de Wiener/GD y algoritmo LMS/SGD¶
Supuestos: Wiener requiere un ambiente estacionario lo cual nos permite calcular \(R_{uu}\) y \(R_{ud}\). En LMS la señal puede ser no estacionaria.
Aprendizaje: En el filtro de Wiener el aprendizaje es determinista. En LMS el aprendizaje viene promediando a nivel de los estimadores de \(w\). En LMS el aprendizaje es estadístico.
Optimalidad: Wiener es óptimo en cambio LMS es sub-óptimo (localmente óptimo). LMS tiende a la solución de Wiener
Costo: LMS se actualiza online y tiene costo \(L\). Wiener se entrena offline y tiene costo \(L^2\)
A continuación se muestra un diagrama que compara el filtro de Wiener y el algoritmo LMS